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Summary. In this work we review some recent development in the mathematical 
modeling of quantitative sociology by means of statistical mechanics. After a short 
pedagogical introduction to static and dynamic properties of many body systems, 
we develop a theory for particle (agents) interactions on random graph. 

Our approach is based on describing a social network as a graph whose nodes 
represent agents and links between two of them stand for a reciprocal interaction. 
Each agent has to choose among a dichotomic option (i.e. agree or disagree) with 
respect to a given matter and he is driven by external influences (as media) and 
peer to peer interactions. These mimic the imitative behavior of the collectivity and 
may possibly be zero if the two nodes are disconnected. 

For this scenario we work out both the dynamics and, given the validity of the 
detailed balance, the corresponding equilibria (statics). Once the two body theory 
is completely explored, we analyze, on the same random graph, a diffusive strategic 
dynamics with pairwise interactions, where detailed balance constraint is relaxed. 
The dynamic encodes some relevant processes which are expected to play a crucial 
role in the approach to equilibrium in social systems, i.e. diffusion of information 
and strategic choices. We observe numerically that such a dynamics reaches a well 
defined steady state that fulfills a shift property: the critical interaction strength 
for the canonical phase transition is lower with respect to the one expected from 
canonical equilibrium. 

Finally, we show how the stationary states of this kind of dynamics can be de- 
scribed by statistical mechanics equilibria of a diluted p-spin model, for a suitable 
non-integer real p > 2. Several implications from a sociological perspective are dis- 
cussed together with some general outlooks. 
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1 Introduction 

Born as a microscopic foundation of thermodynamics, statistical mechanics provides 
nowadays a flexible approach to several scientific problems whose depth and wideness 
increases continuously. In the last decades, in fact, statistical mechanics has invaded 
fields as diverse as spin glasses [MPV87], neural networks [Ami92], protein folding 
[Hua07], immunological memory [Par90], and also made some attempt to describe 
social networks [CG07], theoretical economy [Coo05] and urban planning [CSSOO]. 

In this paper we study statistical mechanics of imitative diluted systems, paying 
particular attention to its applications in social sciences. 

After a review of the statistical mechanics methodology, we introduce the tools, 
both analytical and numerical, used for the investigation of many-body problems. 
We apply such a machinery to study at first the global behavior of a large amount 
of dichotomic agents (i.e. able to answer only 'yes/no" to a given question) whose 
decision making is driven both by a uniform external influence (as the media) and 
by pairwise imitative interactions among agent themselves. In general, the agents 
making up a community are not all contemporarily in contact with each other, 
namely, the network representing the social structure is not fully-connected but 
rather randomly diluted, hence mirroring acquaintances or family relationships. 

Even though refined models as small-worlds graphs have recently been proposed, 
a standard one for such a network is provided by the famous Erdos-Renyi graph. 
The model turns out to be non trivial as, tuning the degree of connectivity and/or 
the strength of interaction, a cooperative state among the agents appears. Moreover, 
we show the existence of a region of the parameter space which is more convenient 
for the global behavior of the society, i.e. it corresponds to a minimum in the free 
energy. 

As social systems do not need to obey Maxwell-Boltzmann distribution, de- 
tailed balance is not strictly required for their evolution, so, after having explored 
the "canonical" 2-body model in full detail, we introduce a more realistic de- 
scription of its free temporal evolution by adopting a diffusive strategic dynamics 
[ABCV06a, ABCV06b]. This dynamic takes into account two crucial aspects, which 
are expected to be effective in the temporal evolution of social systems, i.e. diffu- 
sion of information and strategic choices. We implement it on the same Erdos-Renyi 
graph and study its equilibria. Each agent is selected through a diffusive rule, and 
a flip in its dichotomic status is not weighted "a la Glauber" [Ami92] but rather 
according to a strategic rule which produces the maximum energy gain. We stress 
that this operation involves more-than-two-body effective interactions, as the chosen 
agent interacts both with the first selected one as well as with its nearest neighbors, 
as a whole. This dynamics is shown to relax to a well defined steady state, where 
all the properties of stationarity are recovered [EM90], however the strength of the 
interactions at the critical line is lower, of a few percent, than the expected. The 
whole scenario suggests a "latent" many-body coupling influence, encoded into the 
particular rule for selecting the agents. This is also corroborated by further numer- 
ical analysis. As a consequence, we work out analytically a theory for the randomly 
diluted p-spin model so to fit an effective pgl, which turns out to be p = 2.15, in 
order to match the numerical data available by the dynamics. This result has impli- 
cation both in market trends, as well as in quantitative sociology, where the effective 
interactions always play an important role in decision making [Dur99, McFOl]. 
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The paper is organized as follows: In section 2, for the sake of completeness, 
stochastic dynamics for discrete many body problems is outlined, section 3 deals 
with definitions and introduction to their equilibrium via statistical mechanics. In 
section 4, instead, the model we study is solved in full details with the aim of pre- 
senting both a scenario for these decision makers on random graphs as well as a 
general mathematical method which can be extended by the reader to other models. 
In section 5 our alternative dynamics is introduced and shortly discussed; then fur- 
ther numerical investigations toward a better understanding of a p > 2 behavior are 
presented. 

In section 6 the randomly diluted p-spin model is defined and exploited in all details, 
both analytically (within the cavity field framework) as well as numerically (within 
a Monte Carlo approach). Full agreement is found among the two methods. At the 
end, the last section is left for conclusions: the effective interaction is found and its 
implications analyzed. 

Furthermore, even though the paper is written within a theoretical physics approach, 
remarks concerning the application to quantitative sociology are scattered through- 
out the work. In particular, in the conclusion, a toy application of the outlined 
theory to trades in markets is shown. 



2 A brief introduction to many-body dynamics 

In this section we introduce the fundamental principles of stochastic dynamics used 
to simulate the relaxation to equilibrium of the systems we are interested in. Even 
though for discrete systems two kinds of dynamics are available, parallel and se- 
quential, we are going to deepen only the latter as is the one we will implement 
thought the chapter. Although the topic is well known (see e.g. [LL80, Lig99]), for 
the sake of completeness and in order to offer to the reader a practical approach to 
these models, we present the underlying theory. 

2.1 The model 

Let us consider an ensemble of N agents labelled as i = 1, .., JV. Each agent has two 
possible choices, with respect to a given situation, which are encoded into a variable 
Oi = ±1, say Oi = +1 is "agreement" and viceversa for —1. Each agent experiences 
an external influence (for example by media) which is taken into account by the 
one-body coupling Hi(a; 6) 

N 

ffi(o-;0) = -5>i(t)fc(t), (1) 
i=i 

where #;(£) is the stimulus acting on the i th agent at a given time t and Oi(t) is the 
opinion of the i agent at the same time. 

The interactions among the other agents are encoded into the Ho{o;3) term as 
follows 

N N 

Ho(u; J) = -XX J«(a)(Ti(t)oi-(t). (2) 

i j>i 
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For the moment there is no need to introduce explicitly the dilution of the underlying 
random network as the scheme applies in full generality and may be a benchmark for 
future development by the reader himself. We only stress that Jij(a) is quenched, 
i.e. does not evolve with time, and can be thought of as a symmetric adjacency 
matrix in such a way that a zero entry Jij = means that the agents i and j 
are not in contact with each other, viceversa for Jij = 1 there is a link between 
them. The ratio of connections Ylij Jij/N 2 is tuned by a parameter a, such that 
for ct — > oo the graph recovers the fully connected one, while for a — the graph is 
completely disconnected. Overall, the Hamiltonian defining the model is the sum of 
the two contributes, namely (with a little abuse of notation, thinking at J as (J, 6)) 
H(a; J) — Ho(a; J) + Hi(a; 8). With the signs as they are here, a positive value of 
Jij makes the relevant spins want to line up together, i.e. to share the same opinion, 
and each spin also wants to be aligned with the corresponding external field. 

The investigation of the properties displayed by systems described by this kind 
of Hamiltonian are both analytical and numerical. The former relies on series ex- 
pansions, field theoretical methods, cavity and replica approaches. The latter are 
mainly based on Monte Carlo simulations where we directly simulate the temporal 
evolution of the system in such a way that an expectation value is calculated as 
a time average over the states that the system passes through. However, it must 
be underlined that the Hamiltonian contains no dynamical information, hence we 
have to choose a dynamic for our simulation, namely a rule for changing from one 
state to another during the simulation, which results in each state appearing with 
exactly the probability appropriate to it. Several possibility have been introduced in 
the past, ranging from deterministic, e.g. Q2R dynamics, to stochastic, e.g. Glauber 
algorithm and Wolff algorithm. 



2.2 Transition rates and Markov process 

In statistical mechanics, Maxwell-Boltzmann statistics (hereafter simply "Boltz- 
mann statistics" ) describes the statistical distribution of material particles over var- 
ious energy states in thermal equilibrium, when the temperature is high enough and 
density is low enough to render quantum effects negligible. Now, given the generic 
configuration a — {o"i}i=i,...,jv according to Boltzmann statistics, the value of a 
thermodynamic observable X(/3;J) is given by 

J) = (X(a; J))p = , (3) 

where j3 represents the inverse of the temperature (that sometimes we call "noise" ) , 
i.e. f3 = (fcsT) -1 , being fcs the Boltzmann constant, and the brackets are implicitly 
defined by the r.h.s. of eq.(3). 

Monte Carlo techniques [NB01] work by choosing a subset of states S at random 
from some probability distribution p CT which we specify. Our best estimate of the 
quantity X(a; J) is then given by the so-called estimator {X(a;J)) /3 g: 

(X(a;J)) n * = ^< g6g >^ 1_J . ( 4 ) 
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The estimator has the property that, as the number of sampled states \S\ increases, 
it becomes a more and more accurate estimate of X(/3; J), and, as \S\ — > oo we have 
{X{<r,3)) f>>s ={X(*;J)) fl . 

The choice made for p a is based on the following argument: when in equilibrium, 
the system is not sampling all states in S with equal probability, but according to the 
Boltzmann probability distribution. Hence, the strategy is this: instead of selecting 
the subspace 5* in such a way that every state of the system is as likely to gets 
chosen as every other, we select them so that the probability that a particular state 
s is chosen is p CT = — Z~ 1 e~ ,3H< - c '\ Z being a proper normalization factor by 
now. 

The estimator then simplifies into a simple arithmetic average 

{X(a;J))^^-Lj2 X ^ j y ( 5 ) 

This definition for {X(a; J)) fl g works much better than the one we would obtain 
from a uniform distribution for p a , especially when the system is spending the 
majority of its time in a small number of states. Indeed, the latter will be precisely 
the states sampled most often, and the relative frequency with which we select them 
will correspond to the amount of time the real system would spend in those states. 

Therefore, we need to generate an appropriate random set of states, according 
to the Boltzmann weight p^ q \ In general, Monte Carlo schemes rely on Markov 
processes as the generating engine for the set of states to be used. Let us introduce a 
(normalized) transition probability W[cr, a'] for any pair a, a' of configurations in the 
phase space 5*. Such a set of transition probabilities, together with the specification 
of an initial configuration, allows to construct a Markov chain of configurations, 
S T = (<Ti, <T2, cr T ). The Markov process is chosen in such a way that, when it 
is run for long enough, starting from any state of the system, it will eventually 
produce a succession of states which appear according to the canonical distribution. 
In order to achieve this, two conditions are sufficient: the condition of ergodicity and 
of detailed balance. 

The former is the requirement that it must be possible, for the Markov process, 
to reach any state of the system from any other state, if it is run for long enough. 
Otherwise stated, Vcr, a', 3t : W*[er, a'] is non null, where W*[(r, er'] just represents 
the probability of reaching a' from a vat steps. The ergodic condition is also consis- 
tent with the fact that, in the Boltzmann distribution, every state a appears with 
non-zero probability. On the other hand, notice that this condition does not require 
that W(a,a') 0, V 'a, a '. 

Conversely, the detailed balance condition ensures that, in the limit r — > oo, a 
given configuration a' appears in the Markov chain S T just according to the proba- 
bility distribution pf, . The detailed balance condition requires that the system is 
in equilibrium (the rate of transitions into and out of any state must be equal) and 
that no limit cycles are present. As a result, the detailed balance condition can be 
stated as 

pi^W[a,a']=p^ ) W[a',a}, (6) 

where the l.h.s. represents the overall rate at which transitions from a to a' occur 
in the system, while the r.h.s. is the overall rate for the reverse transition. This 
condition makes the system exhibit time-reversal symmetry at each move and it 
provides a sufficient (but not necessary) condition ensuring that the application of 
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these transition probabilities leads the system to an equilibrium distribution irre- 
spective of the initial state. Now, since we wish the equilibrium distribution to be 
the Boltzmann one, we choose pi- 69 ' = Z^ 1 e~ ,3H (' T '' J \ obtaining 



The constraints introduced so far still leave a good deal of freedom over the 
definition of the transition probabilities. Indeed, the choice of a proper transition 
probability to apply to the system under study is crucial. In fact there is no kinetic 
information in the Hamiltonian given by eq. (1) and by eq. (2), as it only contains 
information about spin orientation and the spatial distribution of lattice sites. It is 
the transition probability which provides a dynamics, i.e. a rule according to which 
the system evolves. 

In the following, we will be especially interested in the so called single-spin- 
flip dynamics, which means that the states involved in the transition only differ 
for the value of a single spin variable. More precisely, in this kind of dynamics, 
given the configuration a, at each time step a single agent i 6 [1,...,N], is ran- 
domly chosen among the N and updated with probability W[a,Fia] to give rise 
to the the configuration FiO, where the N spin-flip operators Fi is defined as 

Fi<7 = Fi{(7l, (7j, CTjv} = {CTI, — <7j, CTjv}. 

The Metropolis and the Glauber dynamics [NB01] are examples of stochastic 
single-spin-flip dynamics. In particular, for the latter one has for the transition 
rates 



where AiH(a; J) = H(F,o; J) - H(a; J). 



3 Equilibrium behavior 

In synthesis, thermodynamics describes all the macroscopic features of the system 
and statistical mechanics allows to obtain such a macroscopic description starting 
from its microscopic foundation, that is, obtaining the global society behavior by 
studying the single agent based dynamics, and then, using Probability Theory, for 
averaging over the ensemble with the weight encoded by the equilibrium distribution 
p^ q \ This scenario is fully derivable when both the internal energy density of the 
system e(/3, a) and the entropy density s((3, a) are explicitly obtained (we are going 
to introduce such quantities hereafter). Then, the two prescription of minimizing the 
energy e (/?,<*) (minimum energy principle) and maximizing entropy s(/3,a) (second 
law of thermodynamics) give the full macroscopic behavior of the system, expressed 
via suitably averages of its microscopic element dynamics. To fulfil this task the 
free energy f(/3,a) = e(/3, a) — /3 _1 s(/3, a) turns out to be useful because, as it is 
straightforward to check, minimizing this quantity corresponds to both maximizing 
entropy and minimizing energy (at the given temperature), furthermore, and this is 
the key bridge, there is a deep relation among it and the equilibrium measure pf , 



W[a,a'] = W[a',a]e 



-P[H(a';l)-H(a;3)] 



(7) 




(S) 



in fact 
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p { : q) ocexpHMTfaJ)), 
f(j3,a)= lim f N (/3,a) = lim -± log V exp(-/3H JV (a; J)), (9) 

N—*oo N-^oo pl\ — ' 

e(y3,a)= lim e N (/3,a)= lim -df3(/3f N (f3,a)), (10) 

N—*-oo N — >oo 

s(/3,a) = hm sjvC3,a)= lim f/jv(/3, a) + ^O^/nW, a))). (11) 

iV— »oo iV— >oo \ / 



So, once explicitly obtained the free energy, equilibrium behavior is solved and the 
whole works once the equilibrium probability distribution is known, provided we use 
the Hamiltonian H(a; J). 

Before proceeding in this derivation, we need some preliminary definitions: 

At first, in the following, it will be convenient to deal with the pressure A((3, a), 

defined as 

A(P,a)= lim A N (j3,a)=-j3 lim f N (P,a); (12) 

N—*oo N—>oc 

we stress that often we are going to consider results in the "thermodynamic limit" 
N — > oo: such procedure allows us to use implicitly several theorem of convergence of 
random variables from Probability Theory and, when the number of agents is large 
enough, the agreement among results at finite iV and results in the thermodynamic 
limit is excellent, as usually the two differ by factors 0(N~ 1 ) or at worse 0(N~ 1 ^ 2 ). 
Let us now further introduce the partition function defined as 

Z N (f3,a)=J2 e ~ PHN{ °- J) =J2P° q) - ( 13 ) 

CT N CT 

As we do not want a sample-dependent theory, using E for the average over the 
quenched variables (i.e. the connectivity), the quenched pressure can be written as 

A N (/3,a) = ^E]nZ N (P,a), 
the Boltzmann state is given by 

w (s(*> J )) = y h ^ Eg( g ' J ) e "' Hw(j;J ' ( 14 ) 

Z N {P, a) ^ 

CT N 

with its replicated form on s replicas defined as 

«(^;J)) = I1^ (s) C9( ct(s) ' J )) ( 15 ) 

s 

and the total average (g) as 

<S> = E[fl(s(ff;J))]. (16) 
Let us introduce further, as order parameters of the theory, the multi-overlaps 

9i~»=^I>F ) -<tf . ( 17 ) 

i — 1 

with a particular attention to the magnetization m = qi = (1/N)'^2 1 ^ =1 a i and to 
the two replica overlap qi2 = (1/N) Yl =i a \ a 1- 
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It is important to stress that the magnetization, which plays the role of the principal 
order parameter (able to recognize the different macroscopic phases displayed by 
the system), accounts for the averaged opinion into the social network, such that if 
(m) = there is no net preference in global decision, while for (m) — ► 1 there is a 
sharp preference toward the "y es " state and viceversa for (m) — * — 1. Analogously (q) 
accounts for similarity among two different "replicas" of the system (two independent 
realization of the adjacency matrix). 

It is easy to check that when [3 —> (or the interaction strength, that is always 
coupled with the noise-), details of the Hamiltonian are unfelt by the agents, which 
will be on average one half up (yes) and one half down (no), giving a null net 
contribution to (m). 

At the contrary when the system is able to experience the rules encoded into the 
Hamiltonian it is easy to see that 

dfN{/3,a) df N {[3,a) 
86— = <^ ' OJ— = <<*<*> 

Averaging over the whole space of choices, we get the macroscopic response of the 
system in terms of the magnetization (m). 

In the thermodynamic limit, further, self-averaging for this order parameter is ex- 
pected to hold, which is expressed via 

lim (m 2 N ) = lim (mjv) 2 , 

JV-.00 N — >oo 

that means that the mean-value of the order parameter is not affected by the details 
of the microscopic structure in the TV — > oo limit (it is an expression of the Central 
Limit Theorem in this framework). 

From a purely thermodynamical viewpoint the equilibrium behavior (the phase 
diagram) is fundamental because it gives both the phase diagram and the critical 
scenario, so to say, the regions in the space of the tunable parameters /3, a where 
the model displays a paramagnetic (independent agent viewpoint) or ferromagnetic 
(collective agent viewpoint) behavior, by which global decision on the whole society 
can not leave aside. 

To obtain a clear picture of the equilibrium of the social network, we use stan- 
dard techniques of statistical mechanics for positive valued interactions, namely the 
smooth cavity field technique [Bar06]. 

For simplicity, as conceptually this does not change the picture, we deal with the 
simpler case Qi = 6 Vi £ [1, TV]. 



4 Equilibrium statistical mechanics of the "2-body" 
model 

The "2-body" model has a long history in physics, having particular importance 
in interaction theories. In fact, from one side, (apart historical problems dealing 
with the deterministic dynamical evolution of the 3-body problem), for a long time 
the interaction in physics were thought of as particle scattering processes and in 
these events the probability of a more-than-2 bodies instantaneous interaction were 
effectively negligible. From the other side, the structure of the 2-body energy is 
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quadratic in its variables and this encodes several information: firstly, from a proba- 
bilistic viewpoint, the Maxwell-Boltzmann probability distribution assumes the form 
of a Gaussian, then, the forces (the derivatives of the energy with respect to its vari- 
able) are linear, such that superposition principle and linear response theory do hold 
and TLC is respected. However, as we will see later, neither of these properties of 
physical systems need to be strictly obeyed in social theories and the interest to 
p > 2 body interactions will arise. 

In this section we consider diluted systems and systematically develop the in- 
terpolating cavity field method [Bar06] and use it to sketch the derivation of a free 
energy expansion: the higher the order of the expansion, the deeper we could go 
beyond the ergodic (agent independent) region. Within this framework we perform 
a detailed analysis of the scaling of magnetization (and susceptibility) at the critical 
line. The critical exponents turn out to be the same expected for a fully-connected 
system. Then, we perform extensive Monte Carlo (MC) simulations for different 
graph sizes and bond concentrations and we compare results with theory. Indeed, 
also numerically, we provide evidence that the universality class of this diluted Ising 
model is independent of dilution. In fact the critical exponents we measured are 
consistent with those pertaining to the Curie- Weiss model, in agreement with ana- 
lytical results. The critical line is also well reproduced. 

The section is organized as follows: after a detailed and technical introduction of the 
model, in Section (4.1) we introduce the cavity field technique, which constitutes 
the framework we are going to use in Section (4.2) to investigate the free energy of 
the system at general values of temperature and dilution. Section (4.3) deals with 
the criticality of the model; there we find the critical line and the critical behavior of 
the main order parameter, i.e. magnetization. Section (4.4) is devoted to numerical 
investigations, especially focused on criticality. 

4.1 Interpolating with the cavity field 

In this section, after a refined introduction of the model, we introduce further the 
cavity field technique on the line of [Bar06]. 

Given TV points and families of i.i.d random variables uniformly dis- 

tributed on these points, the (random) Hamiltonian of the diluted Curie- Weiss model 
is defined on Ising TV-spin configurations a = (ai, . . . , crjv) through 



where is a Poisson random variable with mean £ and a > 1/2 is the connectivity. 
The expectation with respect to all the quenched random variables defined so far 
will be denoted by E, and is given by the composition of the Poissonian average with 
the uniform one performed over the families 




(19) 




1,JV 



(20) 



As they will be useful in our derivation, it is worth stressing the following proper- 
ties of the Poisson distribution: Let us consider a function g : N — > K, and a Poisson 
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variable k with mean aN, whose expectation is denoted by E. 
It is easy to verify that 

E[kg(k)} = aNE[g{k - 1)] (21) 
d aN E\g(k)] =E\g(k + l)-g(k)] (22) 
d 2 (aN) 2E[g(k)] = E[g(k + 2) - 2g(k + 1) + g(k)]. (23) 

Turning to the cavity method, it works by expressing the Hamiltonian of a system 
made of N + 1 spins through the Hamiltonian of N spins by scaling the connectivity 
degree a and neglecting vanishing terms in N as follows 

H N+ i(a) = - ^2 a iv a iv ~ — Vi u Vj„ — y^q-j^q-iv+i (24) 

t^ — l v— 1 I/— 1 

such that we can use the more compact expression 

H??+i(a) ~ Hn(&) + Hn{&)&n+\ (25) 

with 

# w (a) = -^<7i„. (26) 



iV iV-»oo 

lvTT Q ~" Q ' 

^— i 



So we see that we can express the Hamiltonian of N + 1 particles via the one of N 

particles, paying two prices: the first is a rescaling in the connectivity (vanishing in 
the thermodynamic limit), the second being an added term, which will be encoded, 
at the level of the thermodynamics, by a suitably cavity function as follows: let 
us introduce an interpolating parameter t G [0, 1] and the cavity function <?jv(a,t) 
given by 

&(a,/3;t)= lim & N (a,/3;t) (27) 

iV— >-oo 

= lim E ln^ - p = lim E ln ^'frf . 

The three terms appearing in the decomposition (25) give rise to the structure of 
the following theorem which we prove by assuming the existence of the thermody- 
namic limit. (Actually we still do not have a rigorous proof of the existence of the 
thermodynamic limit but we will provide strong numerical evidences in Section 4.4). 

Theorem 1. In the N — > oo limit, the free energy per spin is allowed to assume the 
following representation 

A(a, 0) = In 2 - a dA ^ P) + W{a, (3; t= 1) . (28) 



Proof 

Consider the N + l spin partition function Zjv+i(a, /3) and let us split it as suggested 
by eq. (25) 
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Z N+1 (a, /J) = e- pHN + l(a) ~ e- f3HNl&) -> 3flN{ii) ' TN + 1 (29) 

where the factor two appears because of the sum over the hidden <jjv + i variable. 

Defining a perturbed Boltzmann state u (and its replica product Q = u> X ... X u>) 

as 

= E ^ w}g( ^Xw (<,) . fifow) = n- w (^ (i) )) 

where the tilde takes into account the shift in the connectivity a —> a and multiply- 
ing and dividing the r.h.s. of eq.(29) by Zn{&,(3), we obtain 

Z N+1 (a,/3) = 2Z N (&,f3)Q(e l3 ^=^). (30) 
Taking now the logarithm of both sides of eq.(30), applying the average E and 
subtracting the quantity [ In Zjv+i(a, /?)], we get 

E[ln Z N+1 (a, (3)] - E[ln Z N+1 {&, (3)} = In 2 + E [in f"^® J + <M«, /?; t = 1) 

L Zjv+i(a, p)J 

(31) 

in the large N limit the l.h.s. of eq.(31) becomes 

E[ln Zn+i (a, (3)] - E[ln Z N+ i (a, /?)] = (32) 

( a -Q)A E [ mZjv+1 ( Qi/ 3)] = a— !— j-[ln^+i(a,/3)] = a 9AiV + l(a ' /3) 
att 7V + 1 oa eta 

then by considering the thermodynamic limit the thesis follows. □ 

Hence, we can express the free energy via an energy-like term and the cavity func- 
tion. While it is well known how to deal with the energy-like [ABC08, Gue95, GT04], 
the same can not be stated for the cavity function, and we want to develop its ex- 
pansion via suitably chosen overlap monomials in a spirit close to the stochastic 
stability [AC98, CG05, ParOO], such that, at the end, we will not have the analyt- 
ical solution for the free energy in the whole (a, f3) plane, but we will manage its 
expansion above and immediately below the critical line. To see how the machinery 
works, let us start by giving some definitions and proving some simple theorems: 

Definition 1. We define the t-dependent Boltzmann state uj t as 



Z N , t (a,[3) j-' 



where Z^,t(oi,P) extends the classical partition function in the same spirit of the 
numerator of eq.(33). 



As we will often deal with several overlap monomials let us divide them among two 
big categories: 
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Definition 2. We can split the class of monomials of the order parameters in two 
families: 

• We define "filled" or equivalently "stochastically stable" all the overlap monomi- 
als built by an even number of the same replicas (i.e. q\ 2 , rn 2 , (j , i2Q34q , i234 / )- 

• We define "fillable" or equivalently "saturable" all the overlap monomials which 
are not stochastically stable (i.e. qi2, m, qi2qzi) 

We are going to show three theorems that will play a guiding role for our expansion: 
as this approach has been deeply developed in similar contexts (as fully connected 
Ising model [Bar08a] or fully connected spin glasses [Bar06] or diluted spin glasses 
[GT04], which are the boundary models of this subject) we will not show all the 
details of the proofs, but we sketch them as they are really intuitive. The interested 
reader can deepen this point by looking at the original works. 

Theorem 2. For large N , setting t — 1 we have 

Cj Nt t{oi 1 a i2 ...a irl ) = Wjv+i(ffi 1 ff« 2 ...(7i„o& + i) + O (yj^j (34) 

such that in the thermodynamic limit, if t = 1, the Boltzmann average of a fillable 
multi-overlap monomial turns out to be the Boltzmann average of the corresponding 
filled multi-overlap monomial. 

Theorem 3. Let Qz„ be a fillable monomial of the overlaps (this means that there 
exists a multi-overlap q2n such that q2nQ2n is filled). We have 

lim lim(Q 2 n)t = (q2nQ2n) (35) 

N — >oo t — >1 

(examples: for N -> oo we get (mi) t — ► (mf), (1312)* — > (qia), (912934)4 — > 

(9l2<?34<7l 2 34),) 

Theorem 4. In the N — > 00 limit the averages {■) of the filled monomials are t- 
independent in (3 average. 

For the proofs of these theorems we refer to [ABC08] . 

It is now immediate to understand that the effect of Theorem (2) on a fillable 
overlap monomial is to multiply it by its missing part to be filled (Theorem (3)), 
while it has no effect if the overlap monomial is already filled (Theorem (4)) because 
of the Ising spins (i.e. oj^ +1 = 1 Vn G N). 

Now the plan is as follows: We calculate the f-streaming of the •Z' function in 
order to derive it and then integrate it back once we have been able to express it as 
an expansion in power series of t with stochastically stable overlaps as coefficients. 
At the end we free the perturbed Boltzmann measure by setting t = 1 and in the 
thermodynamic limit we will have the expansion holding with the correct statistical 
mechanics weight. 







Epnw(e 



dt dt 

= 2«3E[ln Q (e PE -i' ^+^0)] _ 2aE[\nCj(e ^ 2J i t 



(36) 



= 2aE[ln u t {e 
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and by the equality e^ CTi o = cosh f3 + <7j sinh /3, we can write the r.h.s. of eq.(36) as 



^d'/' = 25E [ ln ' :) t( cosll / 3 + 'io sinh/3)] = 2alogcosh/3-2oaEpn(l + (Dt(ffi )d)]. 
We can expand the function log(l + U)t6) in powers of 8, obtaining 

9H&,t) nx a ,^(-iT 



dt 



= 2a In cosh /3- 2a V ^ '—8 n (qi,..„ n ) t . (37) 



We learn by looking at eq.(37) that the derivative of the cavity function is built 
by non-stochastically stable overlap monomials, and their averages do depend on 
t making their t-integration non trivial (we stress that all the tillable terms are 
zero when evaluated at t = due to the gauge invariance of the model). We can 
escape this constraint by iterating them again and again (and then integrating them 
back too) because their derivative, systematically, will develop stochastically stable 
terms, which turn out to be independent by the interpolating parameter and their 
integration is straightforwardly polynomial. To this task we introduce the following 

Proposition 1. Let F s be a function of s replicas. Then the following streaming 
equation holds 

= 2a9[Y l (F s at }t,a - s(F s at Q +1 ) t . A ] (38) 

a=l 

a<b a — 1 

1,3 1,3 

+ 2a6» 3 [ ( F 3O-t o- b i0 a- } t . & - s ^(F 3 ct° CT- CT.f +1 )t, a 

a<6<c a<6 

, «(s + l)v^/ j? a s + 1 s+2\ | s ( s + + 2) J „ _ a + l_ a +2_s+3\ i 

H 2\ l^^ sai o a io io H §] ^ s(J *0 *o a »0 K&J 

a — 1 

where we neglected terms O(0 A ). 

For a complete proof of the Proposition we refer to [ABC08] . 
4.2 Free energy analysis 

Now that we exploited the machinery we can start applying it to the free energy. 
Let us at first work out its streaming with respect to the plan (a,/3): 



dA{a,p) _ (H) 1 



di3 



{a} l/=l 



00 TT7 00 / /3(J ' (7 ' \ 

y: *ir(* -1,^1 ^ ) fc ] = o x: »(* - 1, r K /% e «'Y 1 

aE r ( cosh j±fik^jk sinh ff) ) i = aE r u{°ik°jk) + e i 

L w ( cosh /3 + Oi k Oj k sinh /3) J L 1 + w(<7i,.CTj lt )#J 
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by which we get (and with similar calculations for d a A(a,/3) that we omit for the 
sake of simplicity): 

dA(a,/3) _ 



dp 



■a^2(-l) n (l-ef , )0 n - 1 (^,.., n ) (40) 



n=l 



dA(a,f3) 

= In cosh/? - y 9 (gi,..,„) (41) 

da z — ' n 

n — l 

Now remembering Theorem (1) and assuming critical behavior (that we will verify 

a fortiori in sec. (4.3)) we move for a different formulation of the free energy by 
considering the cavity function as the integral of its derivative. In a nutshell the 
idea is as follows: Due to the second order nature of the phase transition for this 
model (i.e. criticality that so far is assumed) we can expand the free energy in terms 
of the whole series of order parameters. Of course it is impossible to manage all 
these infinite overlap correlation functions to get a full solution of the model in the 
whole (a, 0) plane but it is possible to show by rigorous bounds that close to the 
critical line (that we are going to find soon) higher order overlaps scale with higher 
order critical exponents so we are allowed to neglect them close to this line and we 
can investigate deeply criticality, which is the topic of the section. 
To this task let us expand the cavity functions as 

¥(&,l3,t) = J o %dt' (42) 

= 2at log cosh p + f (m) t >, a dt' - \p0 [ (qi2) t ',&dt' + 0{f) 
Jo 1 Jo 

where j3 — 2a9 — » = 2a8 for TV — > oo. Now using the streaming equation as 
dictated by Proposition (1) we can write the overlaps appearing in the expression of 
\P as polynomials of higher order filled overlaps so to obtain a straightforward poly- 
nomial back-integration for the as they no longer will depend on the interpolating 
parameter t thanks to Theorem (1). 

For the sake of simplicity the 5-dependence of the overlaps will be omitted keeping 
in mind that our results are all taken in the thermodynamic limit and so we can 
quietly exchange ct with a in these passages. The first equation we deal with is: 

= /3[(m 2 ) - (m 1 m 2 }t] (43) 

where (mima) is not filled and so we have to go further in the procedure and derive 
it in order to obtain filled monomials: 

d{mim2)t 



dt 

2f3[(m\m 2 )t - (mim 2 m 3 ) t ] + P8[(m 1 m 2 qi2) - 4(mim 2 <?i3)t + S(m 1 m 2 q3i)t 



(44) 



In this expression we stress the presence of the filled overlap (mim2?i2) and of 
(mf rri2}t which can be saturated in just one derivation. Wishing to have an expansion 
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for {m) t up to the third order in 8, it is easy to check that the saturation of the 
other overlaps in the last derivative would carry terms of higher order and so we 
can stop the procedure at the next step 

d{mlm 2 ) t = ^r/ m 2 m 2xi + ^r unfllled terms i + qiq^ ( 45 ) 
dt 

from which integrating back in t 

(m?m 2 )t = f3[{m\m%))t (46) 
inserting now this result in the expression (44) and integrating again in t we find 

(mim 2 }i = f39{m 1 m2qi2)t + f? '' {m\m\)t 2 (47) 
and coming back to (m)t we get 

~" 2 $ ® 2 $ 2 2 3 

(m) t = (3(m }t —{m\m2q\i)t — {m 1 m,2)t (48) 

which is the attempted result. Let us move our attention to {qi 2 )t, analogously we 
can write 

= 2(3[{m iqi2 }t - <m 3 gi2>t] + Milu) - 4{gi 2 gi 3 )t + 3( 9 i 293 4>t] (49) 

and consequently obtain 

{qvz)t = Pe{q 2 12 )t + P 2 (m im 2qi2)t 2 + 0(6» 4 ). (50) 

With the two expansion above, in the N — > oo limit, putting t — 1 and neglecting 
terms of order higher than 6 4 , we have 

&(a,/3,t = 1) = 2alncosh/3+ y(™- 2 } - -^{mlml) ~ —^—{qi 2 ) - —^-{mimzqw) 

(51) 

At this point we have all the ingredients to write down the polynomial expansion 
for the free energy function as stated in the next: 

Proposition 2. A general expansion via stochastically stable terms for the free en- 
ergy of the random two body interacting imitative agent model can be written as 

A(a, (3) = In 2 + a In cosh /3 + y [0 - l) {ml) + (52) 

It is immediate to check that the above expression, in the ergodic region where 
the averages of all the order parameters vanish, reduces to the well known high- 
temperature (or high connectivity) solution [ABC08] (i.e. A(a, j3) = In 2+a log cosh (3). 
Of course we are neglecting 8 6 higher order terms because we are interested in an 
expansion holding close to the critical line, but we are not allowed to truncate the 
series for a general point in the phase space far beyond the ergodic region. 
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4.3 Critical behavior and phase transition 

Now we want to analyze the critical behavior of the model: we find the critical line 
where the ergodicity breaks, we obtain the critical exponent of the magnetization 
and the susceptibility x, which is defined as (x) = (5N[(m 2 ) — (m) 2 ]. 

Let us firstly define the rescaled magnetization £at as £jv = VNttin- By applying 
the gauge transformation <Tj — > Cjcr jv+i in the expression for the quenched average of 
the magnetization (eq. (48)) and multiplying it times TV so to switch to setting 
t = 1 and sending N — > oo we obtain 

<* ?> = W-T) { ^ 2mim2) + 2(01) { ^ 2qi2) + (53) 

by which we see (again remembering criticality and so forgetting higher order terms) 

that the only possible divergence of the (centered and rescaled) fluctuations of the 
magnetization happens at the value = 1 which gives 2a0 = 1 as the critical line, 
in perfect agreement with the expected Landau-like behavior. The same critical 

line can be found more easily by simply looking at the expression (52) as follows: 
remembering that in the ergodic phase the minimum of the free energy corresponds 
to a zero order parameter (i.e.y/ (m?) = 0), this implies the coefficient of the second 
order a(/3') = 4-(/3' — 1) to be positive. Anyway immediately below the critical 
line values of the magnetization different from zero must be allowed (by definition 
otherwise we were not crossing a critical line) and this can be possible if and only 
if a(/3') < 0. Consequently (and using once more the second order nature of the 
transition) on the critical line we must have a{j3') = and this gives again 2a0 — 1. 

Now let us move to the critical exponents. Critical exponents are needed to char- 
acterize singularities of the observables at the critical line and, for us, these indexes 
are the ones related to the magnetization (m) and to the susceptibility (x)- 
We define r = (2«tanh/3 — 1) and we write (m(r)) ~ Go • r s and (x( T )) ~ Go • t 7 , 
where the symbol ~ has the meaning that the term at the second member is the 
dominant but there are corrections of higher order. 

Remembering the expansion of the squared magnetization that we rewrite for com- 
pleteness 

{m2) = W^T) {mlml) + 2(01) < m ™ 12 > + °(0l) ( 54 ) 

and considering that using the same gauge transformation cr; — > (Tjirjv+i on (eq.(50)) 
we have for the two replica overlap the following representation 

(&) = - {f3 f d 2 _ 1) (m 1 m 2qi 2) + 0{6 6 ) (55) 

we can arrive by simple algebraic calculations to write down the free energy, of 
course close to the critical line, depending only by the two parameters (m 2 ) and 

A(a,0) =ln2 + alncosh/?+^(/3'-l)( ? n?>-^- (^-l) <9i2> +0(6 6 ). (56) 

4 48Q \ Zee I 
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Fig. 1. Left panel: Example of Erdos-Renyi random graph made up of 40 nodes and 
with average degree equal to 4. Right panel: Phase diagram: below a c = 0.5 there 
is no giant component in the graph, a c defines the percolation threshold. Above a c 
at left of the critical line the system behaves ergodically, conversely on the right 
ergodicity is broken and the system displays magnetization. 



By a comparison of the formula obtained by deriving A(ot, 0) as expressed by eq.(56) 
and the expression we have previously found (eq. 41), it is immediate to see that we 
have 

^[^(P'-l)(m\)]=6(ml). (57) 
Close to the value f3' = 1, making a change of variable r = j3' — 1 with d a = 28 d r , 
we get 

h if3 ~ 1)<mi) J ~ 2 a7 [r(mi)1 = 2 <mi) + T^T = e{mi) > (58) 

by which we easily obtain 



d(mj) _ dr 
(ml) ~ t 



(m\) ~ r => ^(m'O ~ A = r S (59) 



Therefore we get the critical exponent for the magnetization, 5 = 1/2, which turns 
out to be the same as in the fully connected counterpart [Bar08a], in agreement with 
the disordered extension of this model [Bar06] . 

Again, by simple direct calculations, once we get the critical exponent for the mag- 
netization it is straightforward to show that the susceptibility (\) obeys 

<X>~M- 1 =t'» (60) 

close to the critical line, by which we find its critical exponent to be once again in 
agreement with the classical fully connected counterpart [Ami92], i.e. 7 = — 1. 
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4.4 Numerics 

In this section, by performing extensive Monte Carlo simulations with the Glauber 
algorithm [NB01], we analyze, Ironi the numerical point ol view, the 2-body diluted 
imitative system previously introduced. 

The Erdos-Renyi random graph is constructed by taking N sites and introducing 
a bond between each pair ol sites with probability p = a/(N — 1), in such a way 
that the average coordination number per node is just a. Clearly, when p = 1 the 
complete graph is recovered. 

The simplest version ol the diluted Curie- Weiss Hamiltonian has a Poisson variable 

per bond as Hn = — X/i/=o a iu a iu i and it is the easiest approach when dealing 
with numerics. 

For the analytical investigation we choose a slightly changed version (see eq.(19)): 
each link gets a bond with probability close to a/N for large N; the probabilities ol 
getting two, three bonds scale as 1/N 2 , 1/N 3 therefore negligible in the thermody- 
namic limit. 

Working with directed links (as we do in the analytical framework) the probability 
ol having a bond on any undirected link is twice the probability for directed link 
(i.e. 2a/N). Hence, for large N, each site has average connectivity 2a. Finally in 
this way we allow sell-loop but they add just cr-independent constant to Hn and 
are irrelevant, but we take the advantage ol dealing with an Hn which is the sum 
ol independent identically distributed random variables, that is uselul for analytical 
investigation. 

When comparing with numerics consequently we must keep in mind that a = 2a. 
In the simulation, once the network has been diluted, we place a spin Oi on each node 
i and allow it to interact with its nearest-neighbors. Once the external parameter /? 
is fixed, the system is driven by the single-spin dynamics and it eventually relaxes 
to a stationary state characterized by well-defined properties. More precisely, alter a 
suitable time lapse to and for sufficiently large systems, measurements of a (specific) 
physical observable x(a,a,f3) fluctuate around an average value only depending on 
the external parameters (i~ x and a. 

Moreover, for a system (a, /3) ol a given finite size N the extent of such fluctuations 
scales as with the size of the system. The estimate of the thermodynamic 

observables (a;) is therefore obtained as an average over a suitable number of (un- 
correlated) measurements performed when the system is reasonably close to the 
equilibrium regime. 

The estimate is lurther improved by averaging over different realizations ol the same 
system (a,/3). In summary, 



where a(t) denotes the configuration ol the magnetic system at time step t and T is 
the decorrelation parameter (i.e. the time, in units ol spin flips, needed to decorrelate 
a given magnetic arrangement). In general, statistical errors during a MC run in a 
given sample result to be significantly smaller than those arising from the ensemble 
averaging. 

Figure (2) shows the dependence ol the macroscopic observables (m) and (e) 
from the size of the system; values are obtained starting from a ferromagnetic ar- 
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Fig. 2. Finite size scaling for the magnetization and the internal energy (inset) for 
a — 10 and 4 = 1.67. All the measurements were carried out in the stationary 
regime and the error bars represent the fluctuations about the average values. We 
find good indication of the convergence of the quantities on the size of the system 
and thus of the existence of the thermodynamic limit. 



rangement, at the normalized inverse temperature /3/a = 1.67. Notice that at this 
temperature the system composed of N = 10 4 is already very close to the asymptotic 
regime. Analogous results are found for different systems (a,/3), with j3 far enough 
from f3 c . 

In the following we focus on systems of sufficiently large size so to discard finite 
size effects. For a wide range of temperatures and dilutions, we measure the average 
magnetization (m) and energy (e), as well as the magnetic susceptibility (x) 

Their profiles display the typical behavior expected for a ferromagnet (i.e. im- 
itative coupling) and, consistently with the theory, highlight a phase transition at 
well defined temperatures f5 c (a). 

Now, we investigate in more detail the critical behavior of the system. We collect 
accurate data of magnetization and susceptibility, for different values of a and for 
temperatures approaching the critical one. These data are used to estimate both 
the critical temperature and the critical exponents for the magnetization and sus- 
ceptibility. In Fig. (3) we show data as a function of the reduced temperature 
r = (|/3 — /3 C |//3 C ) _1 for a = 10 and a = 20. The best fit for observables is the 
power law 

(m) ~ t s , (3 > f3 c (61) 
<X> ~ r'. (62) 

We obtain estimates for /3 c (a), 5(a) and 7(a) by means of a fitting procedure. 
Results are gathered in Tab. 1. Within the errors (< 2% for j3 c and < 5% for the 
exponents), estimates for different values of a agree and they are also consistent 
with the analytical results exposed in Sec. (4.3) 
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Fig. 3. Log-log scale plot of magnetization (main figure) and susceptibility (inset) 
versus the reduced temperature r = (\/3 — /3 C \/ fie)' 1 for a = 10. Symbols represents 
data from numerical simulations performed on systems of size N = 36000, while 
lines represent the best fit. 



a 


Pa' 


5 


7 


10 


9.93 


0.48 


-0.97 


20 


19.92 


0.49 


-1.04 


30 


29.98 


0.48 


-1.04 


40 


39.59 


0.50 


-1.02 



Table 1. Estimates for the critical temperature and the critical exponents S and 
7 obtained by a fitting procedure on data from numerical simulations concerning 
Ising systems of size N = 36000 and different dilutions (we stress that analytically 
we get 5 — 0.5 and 7 = —1). Errors on temperatures are < 2%, while for exponents 
are within 5%. 



We also checked the critical line for the ergodicity breaking, again finding optimal 
agreement with the criticality investigated by means of analytical tools. 



5 Beyond detailed balance: Diffusive strategic dynamics 

As we saw in Section 2, a dynamics which obeys detailed balance is constrained 
to reach the Maxwell-Boltzmann equilibrium. In social context however, detailed 
balance actually lacks a clear interpretation, or a data comparison. In this section 
we figure out a possible and more realistic dynamics which aims to mimic opinion 
spreading and strategic choices on random networks, without any a priori constraint 
based on detailed balance. 
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5.1 The algorithm 

In order to simulate the dynamical evolution of Ising-like system as our one several 
different algorithms have been introduced. In particular, a well established one is the 
so-called single- flip algorithm, which makes the system evolve by means of successive 
spin- flips, where we call "flip" on the node j the transformation aj — ► — Oj [Lig99]. 

More precisely, the generic single-flip algorithm is made up of two parts: first 
we need a rule according to which we select a spin to be updated, then we need 
a probability distribution which states how likely the spin-flip is. As for the latter, 
following the Glauber rule, given a configuration a, the probability for the spin-flip 
on the j-th node reads off as 

^ J) = 1 + e /Mtf(^,J)' ^ 

where AH(a,j,J) = 2tJi J^. JijOj is the variation in the cost function due to the 
flip Oj — > — Oj. Hence, for single-flip dynamics the cost variation AH, consequent 
to a flip, only depends on the spin of a few sites, viz. the j-th one undergoing the 
flipping process and its etj nearest-neighbors. 

As for the selection rule according to which sites are extracted, there exist several 
different choices, ranging from purely random to deterministic. In several contexts 
(condensed- matter physics [BBCV02], sociology [ABC09], etc.) unless no peculiar 
mechanisms or strategies are at work, the random updating seems to be the most 
plausible. 

In a social context a spin-flip can occur as a result of a direct interaction (phone 
call, mail exchange, etc.) between two neighbors and if agent i has just undergone 
an opinion-flip he will, in turn, prompt one out of his on neighbors to change opinion 
where, we recall, in social context opinion plays the role of the spin orientation in 
material systems. 

These aspects are neglected by traditional dynamics and can not be described by 
a random updating rule. A different relaxation dynamics, introduced and developed 
in [ABCV05, ABCV06a, ABCV06b], is able to take into account these aspects, 
namely: 

i. the selection rule exhibits a diffusive character: The sequence of sites selected 
for the updating can be thought of as the path of a random walk moving on the 
underlying structure. 

ii. the diffusion is biased: The 014 neighbors are not equally likely to be chosen but, 
amongst the on neighbors, the most likely to be selected is also the most likely to 
undergo a spin-flip, namely the one which minimizes AH(cr,j,J). 

Let us now formalize how the dynamics works. Our MC simulations are made 
up of successive steps [ABC09] : 

- Being i the newest updated spin/agent (at the very first step i is extracted ran- 
domly from the whole set of agents), we consider the corresponding set of nearest- 
neighbors defined as Mi = {ii, i%, we possibly consider also the subset 
A/i C Mi whose elements are nearest- neighbors of i not sharing the same orien- 
tation/opinion: j G M <3> j G A/i A (Ji<jj = — 1. Now, for any j G A/i we compute 
the cost function variation AH(a,j,J), which would result if the flip Oj — > — (Tj 
occurred; notice that AH(cr,j,J) involves not only the nearest-neighbors of i. 
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- We calculate the probability of opinion-flip for all the nodes in Mi, hence 
obtaining p(a, i\, J),p(a, 12, J), ...,p(er, i ai , J), where p(a,a'j,J) (see eq. 63), is the 
probability that the current configuration a changes due to a flip on the j-th site. 

- We calculate the probability V s (a; i, j; J) that among all possible on opinion- 
flips considered just the j-th one is realized; this is obtained by properly normalizing 
the p(a,j, J): 

V s (*;i,r,J)= ^'/' J) ■ (64) 
2^ J ) 

We can possibly restrict the choice just to the set Mi, hence defining V s (er; J) = 

P( cr ,J- J )/EfeeAi' ! P( ,J ' fc ' J )- 

- According to the normalized probability V s (see eq. 64), we extract randomly 
the node j £ Mi and realize the opinion flip <Tj —> —<Tj. 

- We set j = i and we iterate the procedure. 

Finally, it should be underlined that in this dynamics detailed balance is explic- 
itly violated [BBCV02, ABCV05]; indeed, its purpose is not to recover a canonical 
Boltzmann equilibrium but rather to model possible mechanism making the system 
evolve, and ultimately, to describe, at an effective level, the statics reached by a 
"socially plausible" dynamics for opinion spreading [ABC09]. 



5.2 Equilibrium behavior 

The diffusive dynamics was shown to be able to lead the system toward a well defined 
steady state and to recover the expected phase transition, although the critical 
temperature revealed was larger than the expected one [BBCV02]. Such results were 
also shown to be robust with respect to the the spin magnitude [ABCV05] and the 
underlying topology [ABC09]. More precisely, after a suitable time lapse to and for 
sufficiently large systems, measurements of a (specific) physical observable X(a,f3) 
fluctuate around an average value only depending on the external parameter (5 and 
on the geometry of the underlying structures (in particular, for diluted systems, 
on a). It was also verified that, for a system of a given finite size N, the extent 
of such fluctuations scales as N~z (see also [BBCV02, ABCV05]), as indicated by 
standard statistical mechanics for a system in equilibrium. The estimate of the a 
given observable (X{f3)) can therefore be obtained as an average over a suitable 
number of (uncorrelated) measurements performed when the system is reasonably 
close to the equilibrium regime. 

Moreover, the final state obtained with the diffusive dynamics is stable, well 
defined and, in particular, it does not depend on the initial conditions, thus it 
displays all the properties of a stationary state. Similar to what happens with the 
usual dynamics, the relaxation time needed to drive the system sufficiently close to 
the equilibrium situation is found to depend on the temperature. More precisely, we 
experience the so called 'critical slowing down" : the closer T to its critical value, the 
longer the relaxation time. 

In particular, it was evidenced that there exists a critical value of the param- 
eter pf below which the system is spontaneously ordered. However, /3f was found 
to be appreciably smaller than the critical value /3 c (a) expected for the canonical 
Ising model on a Erdos-Renyi random graph. For example, we found j3f « 0.07 and 
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Fig. 4. Critical behavior of the magnetization for the two dynamics (diffusive and 
standard Glauber) as a function of /3 and fixed a — 10. The former dynamics gives 
rise to a critical point higher with respect to the latter. 



Pf ~ 0.016 for a = 10 and a — 45, respectively, versus /3 C («) = tanh _1 (l/a), yield- 
ing j9c(10) « 0.10 and /3 C (45) = 0.022, (see Fig.4). Interestingly, it is not possible 
describe the system subjected to the diffusive dynamics by introducing an effective 
Hamiltonian obtained from eq. (19) by a trivial rescaling. In fact, calling E the 
numerical energies (to separate them from the analytical e), we consider the depen- 
dence on the magnetization displayed by the energies E s (m) and E(m), measured 
for system evolving according to the diffusive dynamics and to a traditional dynam- 
ics, respectively. As for the latter, from eq. (19) it is easy to see that E oc m 2 . As for 
E s (m), we found that E s < E for < m < 1, while E s = E for m = and m = 1. 
This is compatible with a power law behaviour E s ~ m 2+e . In order to obtain an 
estimate for e we measured the ratio E s / E as a function of m; data are shown in 
the log- log scale plot of Fig. (5) and fitting procedures suggest that e w 0.15. In the 
next section we will show that this result can be interpreted as a consequence of an 
effective p-agent interaction, being p > 2. 



6 Statics of many body interactions 

So far, summing the discussion starting the sections 4 and 5, we understood that 
interactions in social networks may involve more than only couple exchanges. Con- 
sequently, corresponding to this observation the need for a many-body imitative- 
behavior Hamiltonian on a random graph appears. In this section, we introduce 
such a p-spin model and study its thermodynamics. 
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Fig. 5. We extrapolate from m(/3) and E((3) the plot E(m) for both the dynamics. 
It is worth noting that the diffusive dynamics deals to a curve living always "below" 
the one obtained by Glauber dynamics results. This strongly suggest the cooperation 
of p > 2 spins each interaction. 



6.1 The diluted p-agent imitative behavior model 

We start exploiting the properties of a diluted even p-spin imitative model: we 
restrict ourselves only to even values of p for mathematical convenience (the in- 
vestigation with the cavities is much simpler), but, due to monotonicity of all the 
observables in p, there is no need to see this as a real restriction (as simulations on 
odd values of p confirm and we will see later on). 

First of all, we define a suitable Hamiltonian acting on a Erdos-Renyi random graph, 
with connectivity a, made up by N agents <7; = ±1, i £ 

Introducing p families {il}, {it}, of i.i.d. random variables uniformly dis- 

tributed on the previous interval, the Hamiltonian is given by the following expres- 
sion 

H N (a,~/(a)) = - ^ a i i(T i 2...cr i p (65) 

v=l 

where, reflecting the underlying network, k is a Poisson distributed random variable 
with mean value 7(a)7V. The relation among the coordination number a and 7 is 
7 oc q p_1 : this will be easily understood a few line later by a normalization argument 
coupled with the high connectivity limit of this mean field model. 

The quenched expectation of the model is given by the composition of the Pois- 
sonian average with the uniform one performed over the families {i v } 

~ e -^ N h(a)N) k 
EH = E P EM = £ e - Y. H- ( 66 ) 

fc=o ij....<g 

The Hamiltonian (66), as written, has the advantage that it is the sum of (a ran- 
dom number of) i.i.d. terms. To see the connection to a more familiar Hamiltonian 
wrote in terms of adjacency tensor A^,...^ , we note that the Poisson-distributed 
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total number of bonds obeys P 7 jv = 7-/V + 0(v JV) for large N. As there are N p 
ordered spin p-plets (i±, ...,i p ), each gets a bond with probability ~ a/N for large 
JV. The probabilities of getting two, three (and so on) bonds scale as l/N 2 , 1/iV 3 , . . . 
so can be neglected. The probability of having a bond between any unordered p-plet 
of spins is p\ as large, i.e. 2a/N for p = 2. 

It is possible to show that our version of the Hamiltonian in fact is thermody- 
namically equivalent with the more familiar involving the explicit adjacency tensor 
Ai lt ... t i , by recall at first the latter model too: 

JV 

—Hm(o", k) ~ —Hm(o') = ^2 Ai u ... >ip a il ...(T ip , (67) 

ii,...,ip 

where k is a Poisson variable with mean 7iV ~ a p ~ 1 N and Ai lt ... t i k are all indepen- 
dent Poisson variables of mean -y/N p ^ 1 ~ (a/JV) p . 

Then, it is enough to consider the streaming of the following interpolating free 
energy (whose structure proves the statement a priori by its thermodynamics mean- 
ing) , depending on the real parameter t £ [0, 1] 

a 

where k is a Poisson random variable with mean yNt and Ai lt .,.,i are random Pois- 
son variables of mean (1 — t)j/N p ~ 1 . Note that the two models alone are recovered in 
the two extremals of the interpolation (for t = 0, 1). By computing the i-derivative, 
we get 

1 d^) = Eln(1 + n{ ^_ _ a . Pg)tanhm (6g) 
1 N 

- JV? X) Ml + «K-«,)tanh09)) = O, 

*i *$> 

where the label in i§ stands for a new spin, born in the derivative, accordingly to 
the Poisson property (22); as the io's are independent of the random site indices in 
the t-dependent SI measure, the equivalence is proved. 

Following a statistical mechanics approach, we know that the macroscopic be- 
havior, versus the connectivity a and the inverse temperature j3, is described by the 
following free energy density 

A{a,(3) = lim A N (a,/3) (69) 

iV— t-oo 



lim ^-Eln Vexp(-/3J : J]v(a,7(a))). 



a 

The normalization constant can be checked by performing the expectation value of 
the cost function: 

E[H] = - 7 7Vm p 
E[H 2 ] - E 2 [H] = 7 2 /V 2 [(<& - m p ) + O(i)] , (70) 
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by which it is easy to see that the model is well defined, in particular it is linearly 
extensive in the volume. Then, in the high connectivity limit each agent interacts 
with all the others (a ~ N) and, in the thermodynamic limit, a — > oo. Now, if p — 2 
the amount of couples in the summation scales as N(N — l)/2 and, with 7 = 2a, a 
linear divergence of a (desired to get a finite ratio a/N for each coupling) provides 
the right scaling; if p = 3 the amount of triples scales as N(N — 1)(N — 2)/3! and, 
with 7 = 3!a 2 , again we find the right connectivity behavior. The generalization to 
every finite p < N is straightforward. 

6.2 Properties of the random diluted p-spin model 

Before starting our free energy analysis, we want to point out also the connection 
among this diluted version and the fully connected counterpart (where each agents 
interact with the whole community, not just a fraction as in the random network) . 
Let us remember that the Hamiltonian of the fully connected p-spin model (FC) can 
be written as [Bar08b] 

Hn°(v) = 2NV-1 ffi 1 <T ia ...<T ip , (71) 

l<i 1 <...<ip<JV 

and let us consider the trial function A(t) defined as follows 

A(t) = -Eln^exp [p J2 °il°H-°il +(1-*) V m 1' (72) 

a v 

which interpolates between the fully connected p-spin model and the diluted one, 
such that for t = only the fully connected survives, while the opposite happens for 
t = 1. Let us work out the derivative with respect to t to obtain 

d t A(t) = (p- 1)q p_1 lncosh(/3) (73) 

n 

by which we see that the correct scaling, in order to recover the proper infinite 
connectivity model, is obtained when a — > 00, 13 -> and = 2(p - l)a p_1 tanh(/3) 
is held constant. 

Remark i. It is worth noting that in social modeling, usually, the role of the tem- 
perature is left, or at least coupled together, to the interaction strength J. As a 
consequence, in order to keep j3' fixed, on different network dilution, the strength 
must be rescaled accordingly to 

J = tanh - ( — ■ 1 , 

^{p-l)^- 1 )' 

while, if present, an external field remains unchanged as it is a one-body term, like 
h^2i unaffected by dilution. 

Remark 2. The dilute p-spin model reduces to the fully connected one, in the infinite 
connectivity limit, uniformly in the size of the system. 
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6.3 Properties of the free energy via the smooth cavity approach 

In this section we want to show some features of the free energy corresponding to 
this model, which is investigated by extending the previous method (the smooth 
cavity approach) to the many body Hamiltonian. 

As the generalization is simple and immediate to be achieved by the reader, we skip 
the proofs in this section. 

The starting point is always the representation theorem 

Theorem 5. In the thermodynamic limit, the quenched pressure of the even p-spin 
diluted ferromagnetic model is given by the following expression 

A{a,f3) = In 2 - -B-jLA(a,0) +9(a,0,t = 1), (74) 
p — 1 da 

where the cavity function ^(t, a, /?) is introduced as 



E In 



(75) 

with limjv^cx) 7 = 7 as for a and a in the two body model (see eq. 26) 

#(7,0,t) = Hm # N (j,(3,t). (76) 

N — >oo 



Thanks to the previous theorem, it is possible to figure out an expression for the 
pressure by studying the properties of the cavity function ^(a, 0) and the connec- 
tivity shift d a A(a, P). 
Let us notice that 

4-M a , P) = ip~ l)a P_2 In cosh /3 - (77) 
da 

- (p - l)a p ~ 2 ^— ^-^"«...,n), 
z — ' n 

n=l 

— 9{&,P,t) = 2<5 p ~ 1 lncosh/3- (78) 
dt 

-2&v- 1 ^t±L e n ( q i : i n ) & , t . 

n=l 

So we can understand the properties of the free energy by analyzing the proper- 
ties of the order parameters: magnetization and overlaps, weighted in their extended 
Boltzmann state u>t- 

Further, as we expect the order parameters being able to describe thermodynamics 
even in the true Boltzmann states u, O, accordingly to the earlier Definition (2), 
we are going to recall that filled order parameters (the ones involving even numbers 
of replicas) are stochastically stable, or in other words, are independent by the t- 
perturbation in the thermodynamic limit, while the others, not filled, become filled, 
again in this limit (such that for them cj t —* <-o in the high iV limit and thermody- 
namics is recovered). 
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Theorem 6. In the thermodynamic limit and setting t — 1 we have 

Wjv,t(ffi 1 cri 2 ...<Tj n ) = wjv+i(<7j 1 <7j 2 ...cTj n cr.N--|-i). (79) 

Theorem 7. Let Q a b be a not-filled monomial of the overlaps (this means that 
labQab is filled). We have 

lim \im{Q ab ) t = (qabQab), (80) 

N—>-oa t— >1 

( examples: 

for N — > oo we get {mi)t — > (mf), (912)4 — > (<?i2),)- 

Theorem 8. In the N — > 00 Krnii, t/ie averages (•) 0/ t/ie /i/Zerf polynomials are 
t-independent in (3 average. 

With the following definition 

/3 = 2(p- l)6?- x (81) 
= 2(p-l)a p - 1 1 ^ T e ^2(^-1)^-^ = ^ 

we show the streaming of replica functions, by which not filled multi-overlaps can 
be expressed via filled ones. 

Proposition 3. Let F s be a function of s replicas. Then the following streaming 
equation holds 



d{F s ) t ,s 
dt 



S 

^(F.mr 1 )^ - ^F.mJ+J}*,*] (82) 



a=l 



+ pe[£{F.<fo\ & -s£(F.q£ L _ 1 ) t , ii 

a<b a—l 



Remark 3. We stress that, at the first two level of approximation presented here, 
the streaming has has the structure of a ^-weighted linear sum of the Curie- Weiss 
streaming (9° term) [Bar08a] and the Sherrington-Kirkpatrick streaming (9 1 term) 
[Bar06], conferring a certain degree of independence by the kind of quenched noise 
(frustration or dilution) to mathematical structures of disordered systems. 

Overall the result we were looking for, a polynomial form of the free energy, 
reads off as 

A(a,/3) =ln2 + a p_1 lncosh/3 + (83) 



2 

T 



^( K P'9{ q l^)-( q { 2 )) + 0{9'). 



Now, several conclusions can be addressed from the expression (83): 
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Remark 4- At first let us note that, by constraining the interaction to be pairwise, 
critical behavior should arise [LL80]. Coherently, we see that for p = 2 we can write 
the free energy expansion as 

A(a,(3) p=2 = ln2 + alncosh(/3) - £(1 -/3')<m 2 ) - ^-{q 2 2 ), 

which coincides with the one of the diluted two-body model (eq. 52) and displays 
criticality at 2a8 = f , where the coefficient of the second order term vanishes, in 
agreement with previous results (sec. 4.3). 

Remark 5. The free energy density of the fully connected p-spin model is [Bar08b] 
A(f3') = ln2 + lncosh(/3m p_1 ) — (/3/2)m p , which coincides with the expansion (83) 
in the limit of a — > oo and /3 — > with [3' = 2(p — l)a p ~ 1 8 held constant. 

Remark 6. It is worth noting that the connectivity no longer plays a linear role in 
contributing to the free energy density, as it does happen for the diluted two body 
models [ABC08, GT04], but, in complete generality as p — 1. This is interesting in 
social networks, where, for high values of coordination number it may be interesting 
developing strategies with more than one exchange [NB01]. 

6.4 Numerics 

We now analyze the system lastly described, from the numerical point of view by 
performing extensive Monte Carlo simulations. Within this approach it is more con- 
venient to use the second Hamiltonian introduced (see eq.(67)): 

JV N 

=-/) Cr n A n,---,ip a i2 a i3--- a ip- ( 84 ) 

i-i J2<*3<.--<*p — 1 

The product between the elements of the adjacency tensor ensures that the p — 1 
spins considered in the second sum are joined by a link with i%. 
The evolution of the magnetic system is realized by means of a single spin flip 
dynamics based on the Metropolis algorithm [NB01]. At each time step a spin is 
randomly extracted and updated whenever its coordination number is larger than 
p — 1. For a large enough (at least above the percolation threshold, as obviously it 
is the case for the results found previously) and p = 3, 4 this condition is generally 
verified. The updating procedure for a spin a; works as follows: Firstly we calculate 
the energy variation Aei due to a possible spin flip, which for p = 3 and p = 4 reads 
respectively 

JV 

Aei = 2ai AijA itk aja k , (85) 
j<k=i 

N 

Aei = 2ai AijAi :k Ai }W aja k a w . (86) 

j <fc<U7— 1 



Now, if Aei < 0, the spin-flip Oi — > — Oi is realized with probability 1, otherwise it 
is realized with probability e~ f3A<i . 
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Fig. 6. Binder cumulants Gl{T) for systems of different size N, as shown in the 
legend, and connectivity a — 50 (left panel) and a = 80 (right panel). 



The case p = 3 has been studied in details and some insight is provided also for 
the case p — 4, while for p = 2 we refer to Sec. (4.4). Our investigations concern two 
main issues: 

- the existence of a phase transition and its nature 

- the existence of a proper scaling for the temperature as the parameter a is tuned. 

As for the first point, we measured the so-called Binder cumulants defined as 
follows: 

G*(T(a)) S l-$^-, (87) 

where (■) jv indicates the average obtained for a system of size N [Bin97]. The study of 
Binder cumulants is particularly useful to locate and catalogue the phase transition. 
In fact, at any given connectivity (above the percolation threshold), in the case of 
continuous phase transitions, Gjv(T) takes a universal positive value at the critical 
point T c , namely all the curves obtained for different system sizes N cross each 
other. On the other hand, for a first-order transition Gn(T) exhibits a minimum at 
Tmin, whose magnitude diverges as N. Moreover, a crossing point at T cross can be 
as well detected when curves pertaining to different sizes N are considered. Now, 
Train and T cross scale as T min — T c oc iV -1 and T cross — T c oc iV~ 2 , respectively. 

In Fig. 6 we show data for Gn(T) obtained for systems of different sizes 
(N = 400, N = 500, and N = 800) but equal connectivity (a = 50 and a = 80, re- 
spectively) as a function of the temperature T. The existence of a minimum is clear 
and it occurs for T ~ 625 and T ~ 1600. Similar results are found also for p = 4 and 
they all highlight the existence of a first-order phase transition at a temperature 
which depends on the connectivity a. 

In order to deepen the role of connectivity in the evolution of the system we 
measure the macroscopic observable (m) and its (normalized) fluctuations (m 2 ) — 
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Fig. 7. Magnetization (main figure) and its normalized fluctuations (inset) for sys- 
tems of different sizes and different dilution as a function of /3 a p_1 . The collapse of 
all the curves provides a strong evidence for the scaling of the temperature. 



(m) 2 , studying their dependence on the temperature (3 and on the dilution a. Data 
for different choices of size and dilution are shown in Figure 7. 

The profile of the magnetization, with an abrupt jump, and the correspondent 
peak found for its fluctuations confirm the existence of a first order phase transi- 
tion at a well defined temperature T c whose value depends on the dilution a. More 
precisely, by properly normalizing the temperature in agreement with analytical re- 
sults, namely (3 = (3 a p_1 we found a very good collapse of all the curves considered. 
Hence, we can confirm that the temperature scales like a p ~ l . 

6.5 Diffusive dynamics revisited 

In the previous section we showed that the diffusive dynamics introduced give rise 
to a non-trivial thermodynamic which can not be explained by, say, a temperature 
rescaling. The crucial point is that such a dynamics intrinsically yields effective in- 
teractions which are more-than-two bodies. This idea is supported by two important 
evidences (see 5.2): a larger "critical" temperature f3f and the anomalous power-law 
behavior E s ~ m 2 ' 15 , which suggests the effective p to be 2.15. Now, while p spans 
from two to infinity the critical temperature raises accordingly, hence we want to 
check that there exist a suitable real value of p that matches the critical tempera- 
ture found numerically and that is compatible with the plots of E s (m) and E{m). 
When p — 2 and close to criticality, the temperature for the phase transition is given 
by Pc = tanh _1 (l/2a p_1 ) = tanh _1 (l/2a). For p = 2.15 this expression becomes 
/3 C ~ tanh~ 1 (l/2a p_1 ) = tanh -1 (l/2a 115 ): The ratio among the two expressions, 
when evaluated for a = 10 gets approximately 1.4, in agreement with data depicted 
in Fig. 4. 
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Before concluding we notice that for p > 2,p £ N, ferromagnetic transitions are 
no longer critical phenomena. At the critical line the magnetization is discontinuous 
and a latent heat does exist. However, if p is thought of as real, for p slightly bigger 
than two, as suggested by our data, the "jump" in the magnetization is expected to 
be small and to approach zero whenever p — > 2. 



7 A simple application to trading in markets 

An appealing application of the whole theory developed concerns trading among 
agents: Suppose we represent a market society only with couple exchanges (p — 2), 
then there are just sellers and buyers and they interact only pairwise. In this case if 
the buyer i has money (at = +1) and the seller j has the product (tjj = +1), or if the 
buyer has no money and the seller has no products (<7; = <jj = — 1), the two merge 
their will and the imitative cost function (19) reaches the minimum. Otherwise, if 
the seller has the product but the buyer has no money (or viceversa), their two 
states are different (±) and the cost function is not minimized. In this scenario, the 
random graph connects on average each agent to a acquaintances and this simply 
increases linearly the possibility that each agent is satisfied. In fact, the higher the 
number of "neighbors" , the larger the possibility of trading. 

When switching to the case p = 3, other strategies are available: for example 
the buyer may not have the money, but he may have a valuable good which can be 
offered to a third agent, who takes it and, in change, gives to the seller the money, 
so that the buyer can obtain his target by using a barter-like approach. In this 
case the two frustrated configurations from the previous sketch can be avoided by 
multiplying by a factor Ok = — 1 given by the third contributor fc, or everything can 
remain the same of course if the latter does not agree {ok = +1). Interestingly, we 
find that in this case (p = 3), the amount of acquaintances one is in touch with 
(strictly speaking, the degree of connectivity a) does not contribute linearly as for 
p = 2, but quadratically: this seems to suggest that if a society deals primarily 
with direct exchanges, no particular effort should be done to connect people, while, 
if barter-like approaches are allowed, then the more connected the society is, the 
larger is the satisfaction reached on average by each agent in his specific goal. The 
above scenario, intuitively, seems to match the contrast among the classical barter- 
like approach of villages, where, thanks to the small amount of citizens, their degree 
of reciprocal knowledge is quite high and the money-mediated one of citizens in big 
metropolis, where a real reciprocal knowledge is missing. 



8 Conclusions and Outlooks 

The idea to apply statistical mechanics methods to social and economical sciences 
has appeared several years ago in the history of science. The main drawback of the 
approach is the lack of a proper measure criteria for the utility function. Although 
the comparison may appear somehow risky one can say that the current approach is 
at the same stage of the pre-thermodynamic epoch when it wasn't clear at all that 
heat was a form of energy and both the first and second principle of thermodynamics 
were still to be identified. This parallel was indeed pointed out already by Poincare 
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in reply to the Walras theory of Economics and Mechanics [Wal09]. Poincare said 
moreover that the lack of ability to measure the utility function is not a severe 
obstacle at a preliminary stage. What instead has proved to be a serious problem 
in the development of mathematical approaches to economics and social sciences 
has been the choice of axioms before phenomenology had been studied properly and 
quantitative data had been extensively analyzed (see [Bou08]). The attempt that are 
made nowadays, and we proposed an instance of statistical mechanics nature, are 
to use mathematical models to mimic microscopic realistic dynamics and reproduce 
to some extent the typical macroscopic observed behavior. Further refinements of 
the models will be necessary once the collection of data will start to provide a fit 
with the free parameters introduced and hopefully new principles and axioms will 
emerge after that stage. 
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